Model-based halftoning

ABSTRACT

A model-based printing method and system for generating halftone output images corresponding to gray-scale-coded input signals. Models for individual two-level (e.g., black on white) printer types allow predicted printer error signals to be generated which can be used to modify the gray-scale coded signals in such manner as to produce binary signals which, when applied to the printer, create halftone images of enhanced quality. In an alternative embodiment binary signals are selected which minimize an error function based on the difference between (i) a predicted perceived image corresponding to the gray scale inputs as filtered by an eye-model filter and (ii) the halftone image resulting from filtering of the binary sequence by a filter modeling the printer followed by the eye-model filter.

This application is a continuation of application Ser. No. 08/046,513,filed on Apr. 12, 1993, which is a Continuation Under Rule 1.62 of Ser.No. 07/659,753 filed Feb. 2, 1991 now abandoned.

FIELD OF THE INVENTION

This invention relates to apparatus and methods for generating halftoneimages. More particularly, the present invention relates to suchapparatus and methods using binary-valued picture elements toapproximate a gray-scale image. Still more particularly, the presentinvention relates to printing of images on uniform-colored (e.g., white)paper using selectively placed constant valued spots.

BACKGROUND OF THE INVENTION

Digital Halftoning, sometimes referred to as "spatial dithering", is theprocess of creating a binary approximation to a sampled gray scaleimage. See for example, R. Ulichney, Digital Halftoning, MIT Press,1987. Sampled gray scale values are typically quantized to have one of adiscrete number of values, e.g., 256 or 1024 values. The basic idea indigital halftoning is to replace these quantized picture elements(pixels) from a region of the gray-scale image having an average valueof x (where 0=white and 1=black) with a binary pattern pattern of 1s and0s. In accordance with one halftoning technique, the fraction ofresulting 1s is approximately x. The binary pattern is then convenientlyused with a display device such as a CRT display or a printer to producethe values for the pixels in the gray-scale halftone image. If the 1sand 0s are supplied to a printer where the 1s are printed as black spotsand the 0s are left as white spaces, and if the spots and spaces aresufficiently close together, the eye averages the black spots and whitespaces to perceive, approximately, gray level x. In so perceiving theimage the eye exhibits a low-pass filtering characteristic. The numberof gray-scale samples (pixels) is advantageously equal to the number ofbits in the binary pattern.

Recent years have witnessed increasing demand for digital storage andtransmission of gray-scale images. In part, this is due to theincreasing use of laser printers having a resolution of, e.g., 300 spots(dots) per inch, to produce halftone approximations to gray-scale imagessuch as photographs, art work, design renderings, magazine layouts, etc.The conventional approach to achieving high quality halftone images isto use a high resolution printer. However, it can be shown that theprinter resolution required for transparent halftoning with prior arttechniques is of the order of 1400 dots/inch. Such printers are oftenslow and expensive.

Many prior art halftoning techniques assume that the black area of aprinted binary pattern is proportional to the fraction of ones in thepattern. This means that the area occupied by each black dot is roughlythe same as the area occupied by each white dot. Thus, the "ideal" shapefor the black spots produced by a printer (in response to 1's) would beT×T squares, where T is the spacing between the centers of possibleprinter spots. However, most practical printers produce approximatelycircular spots. It is clear, therefore, that the radius of the dots mustbe at least T/√2so that an all-ones binary pattern is capable ofblackening a page entirely. This has the unfortunate consequence ofmaking black spots cover portions of adjacent spaces, causing theperceived gray level to be darker than the fraction of ones. Moreover,most printers produce black spots that are larger than the minimal size(this is sometimes called "ink spreading"), which further distorts theperceived gray level. The most commonly used digital halftoningtechniques (for printing) protect against such ink spreading byclustering black spots so the percentage effect on perceived gray levelis reduced. Unfortunately, such clustering constrains the spatialresolution (sharpness of edges) of the perceived images and increasesthe low-frequency artifacts. There is a tradeoff between the number ofperceived gray levels and the visibility of low-frequency artifacts.

Other distortions that can occur in commonly used laser printers, suchas the Hewlett-Packard line of laser printers, include the peculiarcharacteristic that a white line surrounded by several black linesappears brighter than when surrounded by two single lines. These causefurther distortions to the perceived gray levels.

Block replacement is one commonly used halftoning technique used toimprove perceived and gray-scale images. Using this technique, the imageis subdivided into blocks (e.g. 6×6 pixels) and each block is "replaced"by one of a predetermined set of binary patterns (having the samedimensions as the image blocks). Binary patterns corresponding to theentire image are then supplied to a printer or other display device.Typically, the binary patterns in the set have differing numbers ofones, and the pattern whose fraction of ones best matches the gray levelof the image block is selected. This block replacement technique is alsoreferred to as pulse-surface-area modulations. See the Ulichneyreference, supra, at pg. 77.

In another halftoning technique known as screening, the gray scale arrayis compared, pixel by pixel, to an array of thresholds. A black dot isplaced wherever the image gray level is greater than the correspondingthreshold. In the so called random dither variation of this technique,the thresholds are randomly generated. In another variation, ordereddither, the thresholds are periodic. More specifically, the thresholdarray is generated by periodically replicating a matrix (e.g., 6×6) ofthreshold values.

A technique known as error diffusion is used in non-printer halftonedisplay contexts to provide halftoning when ink spreading and otherdistortions common to printers are not present. See, for example, R. W.Floyd and L. Steinberg, "An Adaptive Algorithm for Spatial Grey Scale,"Proc. SID, Vol. 17/2, pp. 75-77, 1976.

Like most of the known halftoning schemes, error diffusion makesimplicit use of the eye model. It shapes the noise, i.e., the differencebetween the gray-scale image and the halftone image, so that it is notvisible by the eye. The error diffusion technique produces noise withmost of the noise energy concentrated in the high frequencies, i.e.,so-called blue noise. Thus, it minimizes the low-frequency artifacts.However, since error diffusion does not make explicit use of the eyemodel, it is not easy to adjust when the eye filter changes, for examplewith printer resolution, or viewer distance. Error diffusionaccomplishes good resolution by spreading the dots. It is thus verysensitive to ink-spreading, in contrast to the clustered dot schemeslike "classical" screening. In the presence of ink spreading, errordiffusion often produces very dark images, therefore limiting itsapplication to cases with no ink-spreading.

Model-based halftoning approaches have been described generally in thecontext of printed images. For example, Anastassiou in the paper, "ErrorDiffusion coding for A/D Conversion,", IEEE Trans. Cir. Sys., Vol.CAS-36, No. 9, pp. 1175-1186, September 1989 proposes a "frequencyweighted squared error criterion" which minimizes the squared errorbetween the eye-filtered binary and the eye-filtered original gray-scaleimage. He considers the problem intractable and suggests an approximateapproach based on neural networks. Moreover, the disclosed techniquesassume perfect printing, i.e., printing without distortion. Allebach, inthe paper "Visual Model-Based Algorithms for Halftoning Images," Proc.SPIE, Vol. 310, Image Quality, pp. 151-158, 1981, proposes a visualmodel to obtain a distortion measure that can be minimized, but providesno complete approach to achieve halftoning.

Roetling and Holladay, in the paper "Tone Reproduction and Screen Designfor Pictorial Electrographic Printing," Journal of Appl. Phot. Eng.,Vol. 15, No. 4, pp. 179-182, 1979, propose an ink-spreading printermodel, of the same general type used in the present invention, but usesit only to modify ordered dither so that it results in a uniform grayscale. Since ordered dither produces a fixed number of apparent graylevels, this technique cannot exploit ink spreading to generate moregray levels.

SUMMARY OF THE INVENTION

The above limitations and distortions of prior art halftoning techniquesare overcome and a technical advance provided by the present invention,as will be described below and in the accompanying drawing.

Rather than trying to resist printer distortions, as in the conventionalapproach, the present invention provides methods and apparatus thatexploit such characteristics, thereby increasing apparent gray-scale andspatial resolution. A key element in such methods is therefore anappropriate printer model. The present invention provides a generalframework for such models and some specific models for laser printers.

As noted above, the error diffusion technique has been used primarilywith CRT displays where distortion such as ink spreading is not present.A halftoning technique in accordance with a first aspect of the presentinvention is an adaptation of error diffusion, for use on printers. Moreparticularly, the model-based approaches of the present inventionincorporate a model for printer distortions in a modified errordiffusion technique to produce a substantial improvement overconventional clustered ordered dither, in both spatial resolution andseverity of low-frequency artifacts.

A second aspect of the present model-based invention, uses least-squaresmodel-based halftoning, and includes both a printer model and a model ofvisual perception. It produces an "optimal" halftoned reproduction byfinding the binary image that causes a combination of printer and visualmodels to match (in the sense of minimizing squared error) the output ofthe visual model in response to the original gray-scale image. Forone-dimensional halftoning (by row or column), this method isconveniently implemented using the Viterbi algorithm. This well knownalgorithm is described, e.g., in Forney, G. D., Jr., "The ViterbiAlgorithm," IEEE Proc., Vol. 61, pp. 268-278. This second approachsuccessfully exploits the printer and visual models to produce more graylevels and better spatial resolution than conventional one-dimensionaltechniques.

Another aspect of the present invention permits gray-scale images to beprinted in halftone form while retaining high fidelity and maximumflexibility. Thus the original image is transmitted to any of a varietyof printer locations using gray-scale image encoders, and is halftonedat the receiver, just before printing. This variation of the presentinvention prints facsimile copies to be produced at a variety ofprinters or other output devices with improved accuracy. Apart fromcoding efficiency (the gray-scale values can be sent by optimum codingtechniques), this approach permits the halftoning to be tuned to theindividual printer. The latter is advantageous because printercharacteristics vary considerably, for example, write-black vs.write-white laser printers. In other words it permits model-basedhalftoning to exploit the characteristics of the specific printer.

When the number of gray-scale pixels is not equal to the number spotlocations on a printer, one of several forms of interpolation istypically used to supply intermediate pixel values.

BRIEF DESCRIPTION OF THE DRAWING

FIG. 1 shows a well-known sensitivity characteristic of the human eye.

FIG. 2 shows a simple eye model based on prior art teachings.

FIG. 3 shows an impulse response for a filter used in modeling humanvisual perception.

FIG. 4 shows a frequency response for the filter of FIG. 3.

FIG. 5 shows a pattern of black spots in accordance with an idealprinter model.

FIG. 6 illustrates an ink spreading phenomenon occurring in certainprinters.

FIG. 7 illustrates geometrically, the meaning of certain parameters usedin defining typical printer models.

FIG. 8 illustrates the parameters shown in FIG. 7 for particular dotpatterns.

FIG. 9 illustrates certain aspects of a one-dimensional ink-spreadingprinter model.

FIG. 10 is a block/flow diagram illustrating the prior art errordiffusion halftoning technique.

FIG. 11 is a block/flow diagram illustrating modifications to the errordiffusion techniques illustrated in FIG. 10.

FIG. 12 shows a block/flow diagram of an embodiment of the presentinvention based on least squares error minimization.

FIG. 13 is a block/flow diagram of a facsimile system incorporating anembodiment of the present invention.

DETAILED DESCRIPTION Models of Visual Perception

To help understand the use of printer models in accordance with thepresent invention, a brief introduction to some aspects of human visualperception will be presented. As mentioned above, halftoning worksbecause the eye perceives a set of closely spaced black and white spotsas a shade of gray. Alternatively, it may be said that the eye acts asif it contained a spatial low pass filter.

Numerous researchers have estimated the spatial frequency sensitivity ofthe eye, often called the modulation transfer function (MTF). Typical ofsuch is the following estimate for predicting the subject quality ofcoded images.

    H(f)=2.6(0.0192+0.114f)exp{-(0.114f).sup.1.1 }             (1)

where f is in cycles/degree. See Mannos, J. L. and D. J. Sakrison, "TheEffects of a Visual Fidelity Critereon on the Encoding of Images," IEEETrans. on Info. Th., Vol. IT-20, no. 4, pp. 525-536, July 1974. This MTFband on the Mannos and Sakrison teachings is plotted in FIG. 1. Asindicated by Eq. (1), the eye is most sensitive to frequencies around 8cycles/degree. Others have variously estimated the peak sensitivity tolie between 3 and 10 cycles/degree. The decrease in sensitivity athigher frequencies is generally ascribed to the optical characteristicsof the eye (e.g. pupil size). FIG. 1 shows that the sensitivity of theeye has dropped 3 db from its peak at about 3 and 16 cycles/degree, 20db at 35 cycles/degree and about 46 db at 60 cycles/degree. The decreasein sensitivity at low frequencies accounts for the "illusion ofsimultaneous contrast" (a region with a certain gray level appearsdarker when surrounded by a lighter gray level than when surrounded by adarker) and for the Mach band effect (when two regions with differentgray levels meet at an edge, the eye perceives a light band on the lightside of the edge and a dark band on the dark side of the edge).

The eye is more sensitive to horizontal or vertical sinusoidal patternsthan to diagonal ones. Specifically, it is lest sensitive to 45 degreesinusoids, with the difference being about 0.6 db at 10 cycles/degreeand about 3 db at 30 cycles/degree. This is not considered to be large,but it is used to good effect in the most commonly used halftoningtechnique for printers as will be described more completely below.

Many models have been proposed that attempt to capture the centralfeatures of human visual perception. For example, see Jain, A. K.,Fundamentals of Digital Image Processing, Prentice Hall, EnglewoodCliffs, N.J. 1989, especially pp. 56-57; Cornsweek, T. N., VisualPerception, Academic Press, New York, N.Y., 1970; and Netravali, A. N.,and B. G. Haskell, Digital Pictures: Representation and Compression,Plenum, New York, N.Y., 1988, especially pp. 292-297. The simplestvisual perception models include just a filter, for example the filterof Eq. (1). Another, and perhaps most commonly cited include amemoryless nonlinearity, as shown in FIG. 2. There, the input image,represented by x, is shown being subjected to the memorylessnonlinearity 201 to produce a modified image, y, before being filteredby filter 202, e.g., that of Eq. (1). The output, z, of filter 202 isthe perceived image. Such nonlinearities account for Weber's law, whichsays that the smallest noticeable change in intensity is proportional tointensity (intensity=1-gray level). Most commonly it is represented as alogarithm or power law (e.g., 1-(1-x)^(1/3)). More complex modelsinclude, for example, a filter before the nonlinearity 201 or a bank offilters in place of 202.

In many cases, practical considerations dictate a finite impulseresponse (FIR) filter for modeling eye characteristics. Indeed, forcertain of the least-squares halftoning techniques described below itproves advantageous to use a one-dimensional discrete-space model of theform

    z.sub.k =M(x.sub.k+m, . . . ,x.sub.k+m),                   (2)

where the x_(k) 's are samples of the image (from one line or onecolumn), the z_(k) 's are the model outputs (upon which cognition isbased), and M(.) is a sliding-window function with 2m+1 arguments (m isa non-negative integer). Such a model can easily incorporate amemoryless nonlinearity and an FIR filter. Typical models that can beused are typically of the form

    z.sub.k =n(x.sub.k *h.sub.k),                              (3)

where n(.) is a memoryless nonlinearity, h _(m), . . . , h_(m) is theimpulse response of an FIR filter and * denotes convolution. Alsoappropriate in some circumstances for the nonlinearity function 201 is

    n(x)=1-(-x).sup.r                                          (4)

for various values of r. For example, others have found r=1/3 to bebest. While it is advantageous to choose m as large as possible, a valueof m=7 with a 15-th order FIR filter that roughly matched (1) forsamples taken at 300 dpi and viewed at 30 inches was found to involve areasonable level of complexity for many applications. Approximations tothe resulting impulse and frequency response are shown in FIGS. 3 and 4,respectively. In FIG. 4, the dotted curve shows the eye MTF of FIG. 1for comparison; f_(s) =1/τ=157.1 cycles/degree. The asymmetry of theimpulse response in FIG. 3 is an artifact of the filter design program.In FIG. 3, τ is equal to 0.0064 degrees.

Printer Models, Generally

This section will introduce a framework for printer models and somespecific models for laser printers. A good model is one that accuratelypredicts the gray levels produced by a printer. While the teachings ofthe present invention may be applied to a wide variety of printer types,it proves especially advantageous to employ "write-black" laser printershaving, for example, 300 dpi resolution. Typical of such printers arevarious ones of the family of laser printers marketed by theHewlett-Packard company, or the Model LZR 1260 by Data Products.

To a first approximation, such printers are capable of producing blackspots (more commonly called dots) on a piece of paper, at any and allsites whose coordinates are typically given in the form (iT, jT), fori=1, . . . , N_(H) and j=1, . . . , N_(w), where T is the horizontal andvertical spacing between dots (typically in inches), N_(H) is the numberof printable lines (rows of dots), N_(w) is the width of a printableline in dots, and (iT, jT) are the coordinates of the jth site from theleft, on the ith line from the top. (These coordinates are consistentwith matrix notation rather than the usual convention for the plane.)The reciprocal of T is generally referred to as the "printer resolution"in dots per inch (dpi). The site with coordinates (iT, jT) will in thefollowing description be called "site (i,j)". The printer is controlledby sending it an N_(h) by N_(w) binary array B=[b_(i),j ], where b_(i),j=1 indicates that a black dot is to be placed at site (i,j) and b_(i),j=0 indicates that the site is to remain white. The latter will bereferred to as a "white" dot.

As illustrated in FIG. 5, black dots produced by an "ideal" printer areblack circles (no shading) with radius 0.707 T. The latter is thesmallest radius such that black circles placed at all sites completelycover the page. The area of such a dot is 1.57 T², i.e., 57% larger thana T×T square. Accordingly, horizontally or vertically (but notdiagonally) neighboring black dots overlap, and white dots are darkenedby neighboring black dots. Specifically, if a white dot has dhorizontally or vertically neighboring (contiguous) black dots, then14.3 d % of it is blackened.

With an actual printer the black dots are not perfectly round, they'renot perfectly black, they are not the ideal size, and they may besomewhat misplaced. Other practical considerations apply to real, ratherthan ideal, printers. For example, a white line surrounded by a pair ofblack lines is not as bright as when surrounded by several black lines.There are many potential causes for such distortions, e.g., inkspreading, spreading of the laser beam, interaction of the laser and thecharge applied to the drum, the movement of toner particles in reactionto charge, the heat finishing, reflections of light within the paper,and so on.

It should always be kept in mind that an individual dot at a site (i,j)may only assume one of two values, typically black or white. However, asa result of phenomena such as those mentioned above, the apparent graylevel produced by the printer in the vicinity of site (i,j) depends in acomplicated way on b_(i),j and neighboring bits. Thus, due to the closespacing of dots and the limited spatial resolution of the eye, theapparent gray level can be modeled as having a constant value p_(i),j inthis vicinity. That is, although the gray level is not actuallyconstant, the eye responds, only to an average gray level over the site.It is this average gray level that p_(i),j represents.

In accordance with one aspect of the present invention, therefore, aprinter model takes the general form

    p.sub.i,j =P(W.sub.i,j) 1≦i≦N.sub.H,1<+j≦N.sub.W(5)

where W_(i),j consists of b_(i),j and the bits in its neighborhood andp_(i),j is the apparent gray level in the vicinity of site (ij). Forconcreteness, it is helpful to visualize the model as producing a graylevel at all points in a page (not just at integer sites). In this sensea continuous parameter model analogous to the discrete parameter modelof Eq. (5) is given by ##EQU1## where u(s,t) denotes the model graylevel at a point s inches from the left and t inches down from the topof a page or other defined surface, and ##EQU2##

In tailoring a model of the above form to a given printer, a main taskis to identify how the function P specifying p_(i),j depends on the bitsin the neighborhood of b_(i),j. Though a variety of phenomena cancontribute to this dependence, it proves advantageous from an analysisand computational viewpoint to limit the dependence of p_(i),j to one inwhich p_(i),j is determined by the values of the binary matrix arrayB=[b_(i),j ] in a fixed window around the site (i,j). In an illustrativeembodiment of the present invention, a 3×3 window centered on site (i,j)is conveniently used, though other windows may be appropriate in otherparticular embodiments. With this typical 3×3 window, the possiblevalues of P can be listed in a table, e.g., with 2⁹ elements.

An Ink-Spreading Model

A common distortion introduced by most printers is, as illustrated inFIG. 6, that their dots are larger than the minimal coveting size, aswould occur, e.g., if "ink spreading" occurred. An illustrative"ink-spreading" printer model that accounts for this phenomenon is##EQU3## where W_(i),j denotes the window surrounding b_(i),j consistingof b_(i),j and its eight neighbors, as indexed below, using compassdirections, ##EQU4## Function f₁ is the number of horizontally andvertically neighboring dots that are black (i.e., the number of ones inthe set {b_(n), b_(e), b_(s), b_(w) }) function f₂ is the number ofdiagonally neighboring dots (i.e., among {b_(nw), b_(ne), b_(se), b_(sw)}) that are black and not adjacent to any horizontally or verticallyneighboring black dot (e.g., in FIG. 6, for the identified site (i,j),b_(nw) =1 and b_(n) =b_(w) =0). Function f₃ is the number of pairs ofneighboring black dots in which one is a horizontal neighbor and theother is a vertical neighbor (e.g., b_(n) =b_(w) =1 would be one suchpair). The constants α, β and γ are the ratios of the areas of therespective shaded regions shown in FIG. 7 to T².

In terms of the ratio ρ of the actual dot radius to the ideal dot radiusT/√2we have ##EQU5## The above assumes 1≦ρ≦√2; i.e., the black dots arelarge enough to cover a T×T square, but not so large that black dotsseparated (horizontally or vertically) by one white dot would overlap.The parameter α, which is the largest of the three factors, representsthe fraction of a horizontally or vertically neighboring site covered bya black dot. For convenience, this model will be referred to as the αink spreading model. It should be noted that the model is not linear inthe input bits, due to the fact that paper saturates at black intensity.For an ideal printer (no ink spreading) ρ=1, the minimum value, andα=0.143, β=0 and γ=0. For ρ=√2, the maximum value, α=0.46, β=0.079 andγ=0.21.

For a typical printer of the class noted above ρ≈1.25. This valueresults in α=0.33, β=0.029 and γ=0.98. FIG. 8 illustrates how the dotpattern in FIG. 6 is modeled with these values. To illustrate one use ofthis model using a 3×3 matrix of surrounding values to predict theeffective gray scale in an area, it is useful to consider the array ofbinary values which includes, for each horizontal string, the repeating6-bit patterns shown in the left column in Table 1. For example, onesuch horizontal string would be 100000100000 . . . 100000. Thishorizontal string is then replicated vertically, i.e., identical ones ofsuch strings occur from the top of the image to the bottom of the image.Table 1 illustrates some interesting statistics relating to such anarray.

                  TABLE 1                                                         ______________________________________                                                             Darkness Predicted by                                               Frequency Ink-Spreading Model                                      Pattern    of 1's    Window 3 (α = 0.33)                                ______________________________________                                        100000     .17       .28                                                      100100     .33       .55                                                      101000     .33       .55                                                      110000     .33       .44                                                      101010     .5        .83                                                      101100     .5        .72                                                      111000     .5        .61                                                      110110     .67       .89                                                      101110     .67       .89                                                      111100     .67       .78                                                      111110     .83       .94                                                      111111     1.0       1.0                                                      ______________________________________                                    

Since the selected patterns appearing in Table 1 are horizontallyperiodic, the gray level of a white dot depends only on the presence orabsence of horizontally neighboring black dots. Specifically, the graylevel of a white dot is α, 2 αor 0, depending on whether there are one,two or no horizontally neighboring black dots. One can see from the graylevels predicted in Table 1 that the ink-spreading model does much toexplain how patterns with the same numbers of ones can have differentgray levels. For example, it predicts the relative gray levels among thepatterns with 3 ones. On the other hand it does not explain why thepattern 110110 produces an image which is darker than the pattern101110, or why 101010 produces an image which is darker than 111110.

One-Dimensional Models

It proves convenient for some purposes to adopt a simplifiedone-dimensional printer model. This is equivalent to a model forprinting one line (or column) or as a model for printing vertically (orhorizontally) invariant images, i.e. those having all horizontal (orvertical) lines the same, as for the patterns of Table 1. With such amodel, the input to the printer is a one-dimensional sequence

where

    p.sub.k =P(W.sub.k),                                       (13)

W_(k) denotes the bits in some neighborhood of b_(k) and P(W_(k)) issome function thereof. A one-dimensional version of the ink-spreadingmodel presented above is ##EQU6## where W_(k) =(b_(k-1), b_(k), b_(k+1))is a window surrounding b_(k) and δ is a parameter between 0 and 1. Asillustrated in FIG. 10, this model reflects those situations in which ablack dot overlaps a fraction δ of the neighboring sites to the left andthe right. Again the model output is not linearly related to input bits.

To identify the parameter δ, it proves convenient to view this model asa projection of the two-dimensional ink spreading model onto onedimension. Accordingly, δ=α=0.33 has been found to be a good value fortypical ones of the class of printers noted above. Further discussion ofone-dimensional models will assume δ=α. Note that for horizontally(vertically) periodic patterns, the one-dimensional model predictsexactly the same gray levels as the two-dimensional model withinterleaved horizontal (vertical) all-zero lines.

The need for one-dimensional ink-spreading models with window sizelarger than 3 will become apparent when considering that in Table 1 the101010 pattern appears about as dark as 110110, even though it only hasthree-fourths as many 1's. A close examination of printer output showsthat the white line in the middle of 11011 appears much larger andbrighter than the white dot in the middle of 01010. Moreover, the whitedot in the middle of 1110111 appears larger and brighter than that inthe middle of 0110110. When requirements so dictate such effects can bebetter captured in a printer model with window size larger than 3 (or3×3 in two dimensions). Thus, while the particular windows andparameters used in the illustrative printer models given above areuseful for predicting perceived gray levels with improved accuracy,particular models may require adaptation as dictated by more completeinformation about (and control of) the underlying physical parameters(e.g., extent of ink spreading), or by more complete understanding ofperceptual considerations.

Error Diffusion Halftoning Technique

Error diffusion halftoning techniques have been described generallyabove. To facilitate a better understanding of improvements achieved bythe present invention, some aspects of this prior art technique will nowbe reviewed.

In the error diffusion halftoning each image pixel is compared to athreshold which depends upon "prior" image pixels, usually above and tothe left. Alternatively viewed, each image pixel is compared to a fixedthreshold, after a correction factor is applied to its original graylevel to account for past errors. Let [x_(i),j ] be a two-dimensionalgray-scale image (after possible interpolation to include the samenumber of dots as the desired binary image, such as e.g., the creationof additional input signals when the input gray-scale image is smallerthan the binary image to be printed), where x_(i),j denotes the pixellocated at the j-th row and the j-th column. It is useful to assume thatthe image has been scanned, and will be processed left to right and topto bottom. Other orderings are, of course, possible as necessary inparticular cases. The binary image [b_(i),j ] produced by errordiffusion is obtained by the following set of equations ##EQU7## Herev_(i),j is the "corrected" value of the gray-scale image. The errore_(i),j at any "instant" (i,j) is defined as the difference between the"corrected" gray-scale image and the binary image. The "past" errors arelow-pass filtered and subtracted from the current image value x_(i),jbefore it is thresholded to obtain the binary value b_(i),j, where[h_(i),j ] is the impulse response of the low-pass filter. Thus errorsare "diffused" over the image.

A diagram of the error diffusion algorithm is shown in FIG. 10. Thethreshold t represented by block 110 in FIG. 10 is fixed at theexemplary value 0.5, the middle of the gray-scale range. Differenceelements are shown as 120 and 125 in FIG. 10. Typically, a page image isscanned left to right and top to bottom i.e., starting at the top leftand finishing at the lower right. The low-pass filter h_(i),jrepresented by block 115 in FIG. 10 has non-symmetric half-planesupport, the two-dimensional equivalent of causality. That is, theeffect of a prior pixel (to the left or above) can be accounted for, buta future pixel, not yet having occurred, does not contribute to anyerror signal. The filter coefficients are positive and their sum isequal to one, thereby assuring stability. Error diffusion halftoningusually requires only one pass through the data.

Various error diffusion filters have been suggested in the literature(see the Ulichney paper, supra). In the following examples a filterproposed by Jarvis, Judice and Ninke in "A Survey of Techniques for theDisplay of Continuous-Tone Pictures on Bilevel Displays," Comp. Graphicsand Image Processing, Vol. 5, pp. 13-40, 1976, will be used. The filteris characterized by Table 2.

                  TABLE 2                                                         ______________________________________                                                         •                                                                                     ##STR1##                                                                          ##STR2##                                    ##STR3##                                                                                 ##STR4##                                                                            ##STR5##                                                                                   ##STR6##                                                                          ##STR7##                                    ##STR8##                                                                                 ##STR9##                                                                            ##STR10##                                                                                  ##STR11##                                                                         ##STR12##                                  ______________________________________                                    

                  TABLE 3                                                         ______________________________________                                        •                                                                                         ##STR13##                                                                           ##STR14##                                             ______________________________________                                    

In the one-dimensional version of error diffusion the illustrativevalues to be used for the filter are shown in Table 3. There is nofundamental difference between the one- and two-dimensional versions oferror diffusion.

Use of Printer Models in Halftoning

Through the use of printer models described above, the present inventionovercomes many of the disadvantages of the prior art halftoningtechniques, including those present in error diffusion halftoning. Ablock/flow diagram reflecting one aspect of a modified error diffusionsystem that compensates for ink spreading is shown in FIG. 11. Only theink spreading contributions of the "present" and "past" pixels are used.Images printed using the system represented in FIG. 11, with printermodels characterized by Eq. (8) or Eq. (14) have the improved apparentgray level and, at the same time, have the sharpness characteristic oferror diffusion. In particular, the performance of this modified errordiffusion system in regions of rapidly changing gray level and in thepresence of printer distortions is very good.

In regions of constant gray level, the modified error diffusionalgorithm of the present invention produces at least as many gray levelsas the prior art "Classic" technique. In common with prior errordiffusion techniques, the model-based modifications in accordance withthe present invention minimize low-frequency artifacts by shaping thenoise, i.e., moving it to the higher frequencies where it is not visibleor moving it to a blue noise range, where it proves very pleasant to theeye. In regions of slowly changing gray level, error diffusion does notsuffer from the false contouring; there is no need to add microdither tothe image.

The system of FIG. 11 differs in overall organization from that of theprior art system shown in FIG. 10 by the inclusion and use of theprinter model 140. Thus, in particular, the output of the thresholdingoperation, i.e., the actual binary pattern sent to the printer(represented by the functional block 135 in FIG. 11 ), is no longer usedto generate the error signal to be fed back to modify the input grayscale values before submitting them to the thresholding step. Rather, amodified version of the binary pattern processed in accordance withprinter model 140, and reflecting the particular characteristics of theprinter, is used as the feedback sequence. This printer model mayadvantageously take the form of Eqs. (8-12) or Eq. (5). As in the priorart, this feedback sequence from difference circuit 145 is low passfiltered using, e.g., Eqs. (15-17), with the coefficients of Table 2 orTable 3 above. It will be understood by those skilled in the art thatdifferent particular filtering coefficients may be used. It should benoted that the use of past error values in filter 150 is accomplished instandard fashion by storing required past signals in memory forming partof the digital filter 150. The modified error diffusion algorithm thatcompensates for dot overlap is shown in FIG. 11. The modified errordiffusion equations are ##EQU8## where (m,n)<(i,j) means (m,n) precedes(i,j) in the scanning order and

    P.sub.m,n.sup.i,j =P(W.sub.m,n.sup.i,j) for (m,n)<(i,j)    (21)

where W_(m),n^(i),j consists of b_(m),n and its neighbors, but here theneighbors b_(k),l have been determined only for (k,l)<(i,j); they areassumed to be zero for (k,l)≧(i,j). Since only the dot-overlapcontributions of the "past" pixels can be used in (18), the "past"errors keep getting updated as more binary values are computed.

Listing 1 is a sample computer program in the well-known C Languagewhich, when executed on a typical general purpose computer, e.g., theSpark Station Model 1+ processor marketed by Sun Microsystems, willperform the processing shown in FIG. 11 and described above. Listing 1assumes that the input values for the sampled gray scale image, I_(i),jhave been stored, in the processor's memory as have the low pass filtervalues and other needed data and programs. Those skilled in the art willadapt the procedures in Listing 1 to particular other computers andlanguages as needed. The output to the printer is, as in all casesdescribed herein, the values for b_(i),j.

Particular applications for the above-described printer model-basedhalftoning technique, and those described below, will use otherimplementing hardware and, where appropriate, software to suit thespecific requirements of the application. For example, in a modificationto a printer, the required processing can be accomplished by amicroprocessor incorporated within the printer. Model information andthe controlling software can conveniently be stored in read only memoryunits (ROMs).

FIG. 13 presents a further illustrative embodiment of the presentinvention for a facsimile system. An original gray scale image isscanned by conventional scanner 5 at a transmitting location. Scanner 5produces gray scale pixels for coding by conventional gray scale coder6. Coded gray scale pixels from coder 6 are communicated via channel 2to gray scale decoder 7 at a receiving location. The number of decodedgray scale pixels is then adjusted by interpolation element 8 andprovided to difference element 130. The balance of FIG. 13 resemblesFIG. 11 modified to show the coupling to a printer 9, and the use ofmemory 132 for storage of past corrected gray scale image values, andmemories 141 and 152 within printer model 140 and low pass filter 150,respectively. Element 151 provides low-pass filtering as referencedabove.

Printer Model Based Least Squares Error Halftoning

An alternative to the modified error diffusion algorithm described abovewill now be presented. This alternative approach is based on thewell-known least squares error criterion. In this alternative approach,it will be assumed that a printer model, a visual perception model andan image are given. The cascade of the printer and visual perceptionmodels will be called the perceptual printing model. The least-squaresapproach to model-based halftoning in accordance with this aspect of thepresent invention then finds the binary array (one bit per image pixel)that causes the perceptual printing model to produce an output that isas close as possible (with respect to squared error) to the response ofthe visual perception model to the original image. Rather than simplyassuming the eye is a low-pass filter that averages adjacent bits (as inconventional ordered dither and error diffusion), this method activelyexploits the visual model. While previous techniques are sometimesrobust to (tolerant of) printer distortions (such as resistance to inkspreading), the present inventive method actively exploits printerdistortions to create the best possible halftoned reproduction. Theresult is more apparent shades of gray and better tracking of edges.Note that, since the eye filter is noncausal, the least-square approachis also noncausal. That is, the decisions at any point in the imagedepend on "future" as well as "past" decisions. In error diffusion thedecisions at any point in the image depend only on the "past". It isthis noncausality of the present least-squares approach that helps giveit the freedom to make sharp transitions and better track edges.

The least-squares halftoning technique of the present invention isconveniently implemented using a method based on the Viterbi algorithm.See e.g., A. J. Viterbi, "Error Bounds for Convolutional Codes and anAsymptotically Optimum Decoding Algorithm," IEEE Trans. Inf. Th., vol.IT-13, pp. 260-269, April 1967, and G. D. Forney, Jr., "The ViterbiAlgorithm," Proc. IEEE, vol. 61, pp. 268-278, March 1973. Because of theViterbi algorithm, only one pass through the data is required for theleast-squares approach. The present least-squares halftoning methodusing the Viterbi algorithm will now be described in the context ofone-dimensional image x=(x₀, . . . ,x_(N-1)).

The overall system to be described is shown in FIG. 12. There, the inputto the least-squares halftoning system in accordance with the presentinvention is, a one dimensional image x=(x₀, . . . ,x_(N-1)) applied oninput 225 in FIG. 12. A printer model 210 shown receiving the binaryarray b_(k) or input 235, in FIG. 12 may e.g., be of the type givenabove in connection with Eq. (13), and have a sliding-window form givenby

    p.sub.k =P(b.sub.k-n, . . . ,b.sub.k+n).                   (22)

This window will cause the printer output to be based on binary inputsextending n bits in either direction from the current bit. Finally, aneye model 220 including a memoryless nonlinearity n(x), shown as 230 inFIG. 12, followed by a finite impulse response filter 240 with impulseresponse h _(m), . . . ,h_(m), will be used.

In the least-squares approach the halftoned approximation b₀, . . .,b_(N-1) is sought that minimizes the squared error ##EQU9## where, asillustrated in FIG. 12

    z.sub.k =y.sub.k *h.sub.k =n(x.sub.k)*h.sub.k              (24)

    w.sub.k =v.sub.k *h.sub.k =n(p.sub.k)*h.sub.k              (25)

    p.sub.k =P(b.sub.k-n, . . . ,b.sub.k+n)                    (26)

Again, * indicates convolution.

The boundary conditions are

    b.sub.k =0 for k<m+n,k>N-m-n-1

    x.sub.k =0 for k<m,k>N-m-1.

These boundary conditions guarantee that the perceived images (theresponse of the printer perceptual model to the bits and the response ofthe perceptual model to the original image) are perfectly white for k<0and k>N-1, and gradually darken to the "true" intensity in a border ofm+ndots.

In formulating the minimization so the Viterbi algorithm may beconveniently applied, the approach of G. Ungerboeck, "AdaptiveMaximum-likelihood Receiver for Carrier-modulated Data-transmissionSystems," IEEE Trans. Commun., vol. COM-22, pp. 624-636, May 1974, A. S.Acampora, "Maximum-likelihood Decoding of Binary Convolutional Codes onBand-limited Satellite Channel", National Telecommunications Conference,1976, and in A. J. Viterbi and J. K. Omura, Principles of DigitalCommunications and Coding, McGraw-Hill, New York, 1979, pp. 272-277] isfollowed generally, as it results in fewer computations.

As a first step, it be shown that ##EQU10## From Eq. (27), the squarederror ε is the sum of ∥z∥²,k which does not depend on the bits b₀, . . .,b_(N-1) to be selected, plus the γk's, each depending on a differentsubset of the bits. In the Viterbi algorithm, the minimization of ε issimplified by introducing the notion of state, which is a set of bitsfrom which a γk can be determined. Since γk is a function of u_(k-2m), .. . u_(k) and since each v_(j) is a function of b_(j-n), . . . b_(j+n),the state at time k may be taken to be

    S.sub.k =(b.sub.k-2m-n+1, . . . ,b.sub.k, . . . ,b.sub.k+n),(29)

i.e., it consists of 2m+2n consecutive bits neighboring b_(k) will beconsidered to be the "present" bit and b_(k+n) to be the "most recent"bit. The state has been defined so that γk can be determined fromS_(k-1) and S_(k), so that S_(k) can be determined from S_(k-1) and themost recent bit b_(k+n), and so that S_(k) contains as few bits aspossible. In essence, the state S_(k-1) summarizes all that one needs toknow to determine γk expect the present bit. It follows from Eqs. (27),(28) and the definition of state that ##EQU11## where μ (.,.) is afunction determined by Eq. (28) and from the boundary condition S_(m-1)=(0, . . . ,0).

Since there is a one-to-one correspondence between sequences of bits b₀,. . . ,b_(N-1) and sequences of states S_(m-1), . . . ,S_(N-m-1), onemay minimize ε by finding the state sequence S_(m-1), . . . ,S_(N-m-1)that minimizes Eq. (30), rather than finding the binary sequence thatminimizes Eq. (27). It is then a straightforward matter to derive thebinary sequence from the state sequence.

The Viterbi algorithm is an efficient way to find the minimizing statesequence. Let S denote the set of all possible states (the set of allbinary sequences of length 2m+2n). For each k in the range m,m+1, . . .,N-m-1 and for each state sεS the Viterbi algorithm finds a statesequence S_(m), . . . ,S_(k-1), s (ending at time k in state S_(k) =s)for which ##EQU12## is minimum among all state sequences ending in s attime k. Let σ_(k) (s) denote the minimizing state sequence, and letε_(k) (s) denote the resulting minimum value. Then the state sequencethat minimizes Eq. (30) (i.e., the desired solution) is (σ_(N-m-1)(s*),s*) where s* is the state for which ε_(N-m-1) (s*) is the smallest.

For each k starting with k=m and each s, the algorithm finds ε_(k) (s)and σ_(k) (s) using the recursion:

    ε.sub.k (s)=min.sub.S.sbsb.k-1 {ε.sub.k-1 (S.sub.k-1)+μ(S.sub.k-1,s)}                            (32)

    σ.sub.k (s)=(s.sub.k (S*.sub.k-1),s)                 (33)

where S_(k-1) * achieves minimum in εk(s) and S_(m-1) =(0, . . . ,0).

In regard to the complexity of the algorithm, for any state s there areprecisely two states that can precede it. Thus the minimization in Eq.(30) involves two computations of μ(.,.), an addition and a binarycomparison. If sufficient memory is available, the function μ may beprecomputed and saved as a matrix. In this case, the complexity of thealgorithm, in operations per dot, is proportional to the number ofstates: N_(s) =2^(2m+2n). Thus, complexity increases exponentially withm and n, but is independent of the size of the image.

There are ways to reduce the number of states (complexity), at the costof some suboptimality, i.e., an increase in ε. The state reductionapproach based on the following observations: The state at time k-1 wasdefined in such a way that it contained all bits needed to determinev_(k-2m), . . . ,v_(k-1), which in turn enter into the third term of Eq.(28), namely, ##EQU13## Ordinarily, some of the last terms of H, say H_(2m), . . . H _(2m+t-1), are so small that the corresponding terms ofthe sum, v_(j) H_(k-j), can be dropped without much effect. In thiscase, the state at time k may be redefined as

    S.sub.k =(b.sub.k-2m-n+t+1, . . . ,b.sub.k, . . . ,b.sub.k+n),(35)

so that now there are only 2^(2m+2n-t) possible states.

It will be seen that, when compared with the prior art techniques, e.g.,those of Anastassiou, the present invention does not assume perfectprinting. The addition of a printer model in accordance with theteachings of the present invention provides a major advance in the art.Also, the above-described mean square error process provides a closedform solution that will be useful in a variety of particularapplications to those skilled in the art. It is deemed that thoseskilled in the art with the detailed algorithmic description and programexample for the error diffusion embodiment will be well able toimplement the above-described Viterbi-based implementation of themean-square-error algorithm for generating the binary array forapplication to a printer such as those in class of laser printersidentified above.

It should be understood that the above-described models, window sizes,filter coefficients and other system and method parameters are merelyillustrative. Other particular printer (and eye) models may proveadvantageous in particular circumstances, as will be apparent to thoseskilled in the art.

While it has been assumed that the printer parameters are fixed andknown in advance before any of the processing described above, no suchlimitation is essential to the present invention. That is, it isadvantageous in some circumstances to adjust the printer model toaccount for changes in printer parameters, e.g., over time. Inparticular, it is possible to sense printer parameters such as dot sizeor shape as pan of the printing process. Alternatively, such sensing canbe accomplished, as required, separately from the actual printingprocess, i.e., off line. After such sensing, the processing canincorporate the new printer parameters in all future halftoningoperations. ##SPC1##

We claim:
 1. A method for generating an array of output binary signalssuitable for application to a display device to generate a halftoneimage in response to an array of input signals characterizing a grayscale image, said method comprising the steps of:in response to appliedbinary signals, forming past signals predictive of regions of halftoneimages formed by said display device, said past signal formed based on amodel of said display device, modifying each of a plurality of saidinput signals in response to one or more past error signals, said pasterror signals reflecting differences between past modified input signalsand said past signals predictive of regions of halftone images, andforming a binary signal in response to each of a plurality of saidmodified input signals.
 2. The method of claim 1 wherein the step offorming a binary signal comprises assigning one value to said binarysignal whenever a modified input signal exceeds a threshold value, andassigning another value to said binary signal whenever said modifiedinput signal fails to exceed said threshold value.
 3. The method ofclaim 1 wherein the number of input signals in said array of inputsignals is adjusted to be equal to the number of signals in said arrayof output binary signals.
 4. The method of claim 3, wherein for the caseof the number of input signals in said array of input signals being lessthan the number of signals in said array of output binary signals, thenumber of input signals is adjusted by creating additional inputsignals, the value of said additional input signals being determined byinterpolation based on selected neighboring input signals.
 5. The methodof claim 1 wherein said array of input signals is ordered in accordancewith the respective position in said image, and said binary signals usedto generate said predicted halftone image comprise binary signalscorresponding to one or more prior positions in said order.
 6. A systemfor printing halftone images on a printing surface in response to anarray of input signals characterizing a gray scale image,comprising:means for generating spots at selected ones of regularlyspaced positions on said printing surface in response to an orderedsequence of binary signals, means for modifying each of a plurality ofinput signals in response to one or more past error signals, means forapplying a binary signal to said means for generating in response to amodified input signal, means, responsive to applied binary signals, forforming past signals predictive of regions of halftone images formed bysaid means for generating, said past signals formed based on a model ofsaid means for generating, and means for forming past error signalsreflecting differences between past modified input signals and said pastsignals predictive of regions of halftone images.
 7. The system of claim6 wherein said means for applying assigns one value to said binarysignal whenever said modified input signal exceeds a threshold value andanother value to said binary signal whenever said modified input signalfails to exceed said threshold value.
 8. The system of claim 6 furthercomprising means for adjusting the number of input signals in said arrayof input signals to be equal to the number of said regularly spacedpositions.
 9. The system of claim 6, wherein for the case of the numberof input signals in said array of input signals being less than thenumber of said regularly spaced positions, the means for adjusting thenumber of input signals creating additional input signals, the value ofsaid additional input signals being determined by interpolation based onselected neighboring input signals.
 10. The system of claim 6 whereinthe input signals of said array of input signals are ordered incorrespondence with respective positions in said image, and wherein saidpast signals predictive of regions of said halftone images are based onbinary signals corresponding to one or more prior positions in saidorder.
 11. A facsimile system for printing halftone images on a printingsurface at a second location corresponding to a gray scale image at afirst location, the system comprising:means for receiving at said secondlocation an ordered sequence of gray-scale coded input signalsrepresenting said gray scale image, means for generating spots atselected ones of regularly spaced positions on a printing surface inresponse to an ordered sequence of binary signals, means for modifyingeach of a plurality of input signals in response to one or more pasterror signals, means for applying a binary signal to said means forgenerating in response to a modified input signal, means, responsive toapplied binary signals, for forming past signals predictive of regionsof halftone images formed by said means for generating, said pastsignals formed based on a model of said means for generating, and meansfor forming past error signals reflecting differences between pastmodified input signals and said past signals predictive of regions ofsaid halftone images.
 12. The system of claim 11 further comprisingmeans at said first location for generating said ordered sequence ofgray-scale coded signals.
 13. The system of claim 12 wherein said meansfor generating said ordered sequence of gray-scale coded signalscomprises means for scanning said gray-scale image to form a sequence ofvalues corresponding to sequential locations on said image and means forcoding each of said sequence of values.
 14. A method for communicatingan image for printing comprising the steps of:(a) encoding the image andtransmitting its encoded representation; (b) receiving the encodedrepresentation of the image and decoding it; (c) determining a halftoneimage based on the decoded representation of the image, wherein therepresentation of the image comprises one or more input signals, andwherein the step of determining a halftone image includes the stepsof(1) in response to applied binary signals, forming past signalspredictive of regions of halftone images formed by a printing device inresponse to applied binary signals, said past signals formed based on amodel of said printing device, (2) modifying each of a plurality of theinput signals in response to one or more past error signals, the pasterror signals reflecting differences between the past modified inputsignals and said past signals predictive of regions of halftone images,and (3) forming an output signal comprising a binary signal in responseto each of a plurality of the modified input signals; and (d) printingthe halftone image with use of the one or more output signals.
 15. Themethod of claim 14 wherein said forming of a binary signal comprisesassigning one value to said binary signal whenever said modified inputsignal exceeds a threshold value, and assigning the other value to saidbinary signal whenever said modified input signal fails to exceed saidthreshold value.
 16. The method of claim 14 wherein the number of inputsignals is adjusted to be equal to the number of output signals.
 17. Themethod of claim 14 wherein said input signals are ordered in accordancewith the respective position in said image, and said binary signals usedto generate said predicted halftone image comprise binary signalscorresponding to one or more prior positions in said order.
 18. Themethod of claim 17 wherein said binary signals used to generate saidpredicted halftone image comprises binary signals corresponding to oneor more subsequent positions in said regions of halftone images.